\documentclass[12pt,a4paper]{article}
\usepackage{dsfont}
\usepackage{amssymb}
\usepackage[utf8]{inputenc}
\usepackage[russian]{babel}
\usepackage{amsmath}
\usepackage{graphics} 
\usepackage{float}
\usepackage{color} %% это для отображения цвета в коде
\usepackage{listings} %% собственно, это и есть пакет listings
\usepackage{caption}
\DeclareCaptionFont{white}{\color{white}} 
\DeclareCaptionFormat{listing}{\colorbox{gray}{\parbox{\textwidth}{#1#2#3}}}
\captionsetup[lstlisting]{format=listing,labelfont=white,textfont=white}
\begin{document}

\tableofcontents

\newpage
\section*{Введение}
Широко известно, что фрактальные множества являются математическим обобщением и идеализацией многих реальных объектов. Иначе говоря, многие реальные объекты 
могут быть с некоторой погрешностью приближены фрактальными множествами. Естественно эти множества весьма сложны, у них неоднородная структура, нецелая размерность, 
а также они используют много других идеализированных математических объектов (например, круг, прямой угол и т.п.). 
Благодаря мощной вычислительной технике и развитию науки о вычислениях множество результатов, касающихся фрактальных множеств, было получено с помощью компьютера. 
В данной работе описывается метод, представленный Фояшем и Темамом в \cite{1988PhyD...32..163F} в 1988г. Было предложено аппроксимировать аттрактор с помощью 
алгебраического множества, т.е. множества, являющегося корнями какого-то полинома. Доказано, что для данного аттрактора $\mathbb{A}$ и любого положительного 
$\epsilon$ существует многочлен, который можно выписать в явном виде, значения которого на аттракторе не превосходят по модулю $\epsilon$. Нули этого многочлена можно принять 
за аппроксимацию аттрактора.
\newline
Удобство такого подхода состоит в том, что мы можем выписывать такие многочлены с любой точностью. Кроме того, такой подход к аппроксимации аттракторов и фрактальных 
множеств является весьма гибким и выгодно отличается от классических методов аппроксимации аттракторов гладкими многообразиями тем, что поверхности, порожденные корнями 
многочлена не обязаны быть гладкими или иметь один порядок гладкости во всех своих точках, в отличие от многообразия.
\newline
В данной работе изложен общий алгоритм построения таких аппроксимаций, а также частные случае на примере аттрактора Лоренца и инвариантного множества для 
системы Ресслера.

\section{Общие определения}
Будем рассматривать дифференциальные уравнения вида
\begin{displaymath}
u'(t) = F(u(t)),
\end{displaymath}
где $F$ -- гладкая функция, действующая из $\mathds{R}^d$ в $\mathds{R}^d$, а $u(t) = (u_1(t), \cdots, u_n(t))$ -- функция из $\mathds{R}$ (или какого-то 
интервала на вещественной прямой) в $\mathds{R}^d$. Естественным образом $u' = du/dt$.
\newline
Полагается, что для любого $u_0 \in \mathds{R}^d$ задана задача решения дифференциального уравнения с начальными данными
\begin{displaymath}
u(0)=u_0,
\end{displaymath}
определенная на всей положительной полуоси. Обозначим за $S(t)$ отображение
\begin{displaymath}
S(t): u_0 \in \mathds{R}^d \rightarrow u(t)\in \mathds{R}^d.
\end{displaymath}
Нас интересуют диссипативные системы. В случае таких систем траектории не блуждают случайным образом в пространстве, а со временем притягиваются к какому-то 
поглощающему множеству. Дадим ему строгое определение. 
\newline
Ограниченное множество $B_0\subset \mathds{R}^d$ называется поглощающим множеством, если 
для любого ограниченного множества $B\subset \mathds{R}^d$ существует время $t=t_0(B)$ такое, что 
\begin{displaymath}
S(t)B\subset B_0, \forall t\ge t_0(B).
\end{displaymath}

Введем также понятие инвариантного множества: $B \subset \mathds{R}^d$ называется инвариантным, если
\begin{displaymath}
S(t)B = B\:\forall t\ge 0.
\end{displaymath}

Так же классическим свойством диссипативных систем является наличием глобального компактного аттрактора $\mathds{A}$. Аттрактор является подмножеством поглощающего множества 
и является инвариантным замкнутым поглощающим множеством. Структура аттракторов зачастую очень сложна, часто оказывается, что никаких других способов их 
изучения, кроме как компьютерных симуляций траекторий или определения фрактальных размерностей, нет.

\section{Комплексификация решения}
Для дальнейшего изучения предложенной системы проведем ее комплексификацию, заключающуюся в расширении ее области определения на небольшую полоску комплексной плоскости. 
Ради простоты ограничим класс рассматриваемых уравнений следующим:
\begin{displaymath}
u' + Au + R(u) = 0,
\end{displaymath}
где $A$ -- симметричная положительно определенная матрица, а
\begin{displaymath}
R(u) = R_0 + R_1(u) + R_2(u, u),
\end{displaymath}
где $R_0$ -- вектор в пространстве $\mathds{R}^d$, $R_1$ -- линейный оператор на $\mathds{R}^d$, а $R_2$ -- билинейный. Многие классические диссипативные 
дифферениальные уравнения могут быть представлены в таком виде, в том числе и уравнения Лоренца и Ресслера. 
\newline
Для дальнейшего расширения класса рассматриваемых уравнений можно рассматривать $R$ как сумму $k-$линейных операторов над $\mathds{R}^d$, в данной работе ограничимся 
двумя операторами.
\newline
Непрерывность $R_1$ и $R_2$ обеспечивается существованием следующих констант $c_1$ и $c_2$:
\begin{displaymath}
|R_1(u)|\le c_1|u| \:\forall u\in\mathds{R}^d, \: |R_2(u, v)| \le c_2|u|\cdot|v| \:\forall u,v \in \mathds{R}^d.
\end{displaymath}

Проведем теперь процедуру комплексификации. Так как функция $R$ является аналитической функцией переменной $u$, изначальное уравнение с начальными условиями $u(0)=u_0$ обладает 
решением $u(\zeta)$, которое определено в окрестности вещественной полуоси на комплексной плоскости, а именно в множестве $\Omega$, заданном соотношением
\begin{displaymath}
max\{-\delta_0, -Re \zeta\}\le Im \zeta\le min\{\delta_0, Re \zeta\}
\end{displaymath}
Можно показать, что это решение аналитично в данной области (для какого-то $\delta_0$) и ограничено:
\begin{displaymath}
|u(\zeta)|\le M_0, \forall \zeta \in \Omega
\end{displaymath}
Отсюда имеем
\begin{displaymath}
|Ru(\zeta)|\le K_0 =|R_0|+c_1M_0 + c_2M_0^2, \forall \zeta \in \Omega
\end{displaymath}
Рассмотрим теперь какую-нибудь точку на аттракторе, $u_*\in\mathds{A}$. Известно, что $u_*$ принадлежит некоторой орбите $u=u(t), \forall t\in\mathds{R}$, поэтому 
не умаляя общности можно заключить:
\begin{displaymath}
u_*=u(0)
\end{displaymath}
Как показано выше, функция $t \rightarrow u(t)$ может быть расширена на часть комплексной плоскости по аналитичности. Интегрируя исходное уравнение, имеем:
\begin{displaymath}
u(0) = e^{tA}u(t)-\int_t^0 e^{\tau A}R(u(\tau))d\tau
\end{displaymath}
Устремляя $t\rightarrow -\infty$ и учитывая ограниченность $u$, получаем
\begin{displaymath}
u(0) = - \int_{-\infty}^0 e^{\tau A}R(u(\tau))d\tau.
\end{displaymath}
\section{Процедура аппроксимации}
Обозначим за $\Lambda$ какое-нибудь собственное значение матрицы $A$. В силу положительной определенности матрицы все ее собственные числа вещественны. 
Рассмотрим $Q$ - проектор $\mathds{R}^d$ на подпространство, натянутое на собственные вектора матрицы $A$, отвечающие собственным числам, превышающим $\Lambda$. 
Соответственно за $P$ обозначим проектор на оставшуюся часть пространства, т.е. $P=I-Q$.
\newline
Обозначим
\begin{displaymath}
p = Pu, \:\: q = Qu.
\end{displaymath}
Спроецировав начальное уравнение с помощью проектора $Q$, получим:
\begin{displaymath}
\frac{dq}{dt} = Aq + QR(u)=0
\end{displaymath}
Здесь использована коммутативность проекторов с линейными операторами.
\newline
Усилим положительность матрицы $A$ коэффициентом $k$, который в будущем сможем варьировать.
\begin{displaymath}
\frac{dq}{dt} = (A+kI)q -kq + QR(u)=0
\end{displaymath}
А теперь, проинтегрировав полученное уравнение по вещественной оси времени, получим аналог последней формулы из предыдущего раздела:
\begin{displaymath}
q(0) = -\int_{-\infty}^0 e^{\tau(A+kI)Q}(QR(u(\tau)) - kq(\tau))d\tau.
\end{displaymath}
Теперь наша задача состоит в том, чтобы приблизить правую часть этого равенства какой-либо простой алгебраической функцией, так мы сможем аппроксимировать конкретную 
точку на аттракторе вместе с ее окрестностью.
\newline
Для этого разобьем интеграл на два других интеграла, один по промежутку $(-\infty, -\delta)$, а другой по оставшемуся $(-\delta, 0)$. Рассмотрим и оценим первый из них.
\begin{displaymath}
\left|-\int_{-\infty}^{-\delta} e^{\tau(A+kI)Q}(QR(u(\tau)) - kq(\tau))d\tau\right| \le \left|\int_{-\infty}^{-\delta}e^{\tau(\Lambda+k)}(K_0+kM_0)d\tau\right|
\le
\end{displaymath}
\begin{displaymath}
\le (\Lambda+k)^{-1}(K_0 +kM_0)e^{-\delta(\Lambda+k)}
\end{displaymath}
Здесь использована операторная оценка $e^{\tau(A+kI)Q}$, а так же условия на ограниченность функции $R$ на комплексном участке $\Omega$. 
\newline
Заметим, что теперь интеграл оценен сверху значением, которое мы можем сделать сколь угодно маленьким, выбирая соответствующим образом $k$. А значит, можно считать
\begin{displaymath}
q(0)\approx -\int_{-\delta}^0 e^{\tau(A+kI)Q}(QR(u(\tau)) - kq(\tau))d\tau
\end{displaymath}
Продолжим оценивать этот интеграл. Вспомним теперь об аналитичности функции $u$ и разложим ее в ряд:
\begin{displaymath}
u(\tau) = \sum_{n=0}^{\infty} \frac{\tau^n}{n!}u^{(n)}(0),
\end{displaymath}
разложение верно в окрестности $|\tau| < \delta_0$. Соответственно нам надо брать $\delta < \delta_0$. Рассмотрев теперь частичные суммы
\begin{displaymath}
u_N(\tau)=\sum_{n=0}^{N}\frac{\tau^n}{n!}u^{(n)}(0),
\end{displaymath}
можно легко сделать вывод о том, что величины $u_n = u^{(n)}(0)$ являются алгебраическими выражениями от $u_0=u(0)$. Действительно, дифференцируя исходное уравнение и зная предыдущее 
соотношение, можно построить алгебраическую функцию. связывающую $u_0$ и $u_n$.
\newline
Аналогично можно поступить с величинами $f_n = f^{(n)}(0) = R(u(\tau))^{(n)}(0)$, а в силу того, что $f_0$ линейно выражается через $u_0$ с помощью исходного уравнения, 
получаем, что $f_n$ так же являются алгебраическими выражениями от $u_0$.
\newline
Таким образом мы теперь научились приближать функции $u(\tau)$ и $f(\tau)$ алгебраическими величинами. Запишем это следующим образом:
\begin{displaymath}
-\int_{-\delta}^0 e^{\tau(A+kI)Q}(QR(u(\tau)) - kq(\tau))d\tau = -\int_{-\delta}^0 e^{\tau(A+kI)Q}Q(f(\tau) - ku(\tau))d\tau = 
\end{displaymath}
\begin{displaymath}
= -\int_{-\delta}^0 e^{\tau(A+kI)Q}Q(f_N(\tau) - ku_N(\tau))d\tau -\int_{-\delta}^0 e^{\tau(A+kI)Q}Q(f(\tau) - f_N(\tau) + ku_N(\tau) - ku(\tau))d\tau
\end{displaymath}
Второй интеграл можно сделать сколь угодно маленьким, выбрав подходящие значения для $N$ и $k$.
\newline
Запишем теперь первый интеграл в следующем виде:
\begin{displaymath}
-\int_{-\delta}^0 e^{\tau(A+kI)Q}Q(f_N(\tau) - ku_N(\tau))d\tau = \sum_{n=0}^{N}S_nQ(f^{(n)}(0) - ku^{(n)}(0)),
\end{displaymath}
где
\begin{displaymath}
S_n = -\int_{-\delta}^{0}e^{\tau(A+kI)Q}\frac{\tau^n}{n!}d\tau
\end{displaymath}
Проинтегрировав по частям, получаем:
\begin{displaymath}
S_n = (-1)^{n+1}((A+k)Q)^{-n-1} + e^{-\delta(A+kI)Q}\sum_{j=0}^{n}(-1)^{j}\frac{\delta^{n-j}}{(n-j)!}((A+kI)Q)^{-j-1}
\end{displaymath}
Заметим, что первое из этих слагаемых легко можно сделать маленьким с помощью выбора $k$, а второе легко вычислить с помощью компьютера.
\newline
Рассмотрим теперь $u_n$ и $f_n$, для них нам надо предоставить полиномиальную зависимость от $u_0$.
Дифференцируя исходное уравнение, можно записать:
\begin{displaymath}
u_1 = -Au_0 - R(u_0)
\end{displaymath}
\begin{displaymath}
u_2 = -Au_1 - R_1(u_1) - R_2(u_1, u_0) - R_2(u_0, u_1)
\end{displaymath}
\begin{displaymath}
u_3 = -Au_2 - R_1(u_2) - R_2(u_2, u_0) - 2R_2(u_1, u_1) - R2(u_0, u_2)
\end{displaymath}
\begin{displaymath}
\cdots
\end{displaymath}
Далее зависимость очевидна: в $R_1$ необходимо передать предыдущее вычисленное значение, а коэффициенты при $R_2$ распределены по биномиальному закону.
Аналогично для $f_n$ имеем
\begin{displaymath}
f_0 = R(u_0)
\end{displaymath}
\begin{displaymath}
f_1 = R_1(u_1) + R_2(u_1, u_0) + R_2(u_0, u_1)
\end{displaymath}
\begin{displaymath}
f_2 = R_1(u_2) + R_2(u_2, u_0) + 2R_2(u_1, u_1) + R2(u_0, u_2)
\end{displaymath}
\begin{displaymath}
\cdots
\end{displaymath}
Виды самих многочленов можно получить с помощью символьных вычислений, однако сами записи получаются очень громоздкими. 
\newline
Ну и наконец можно разложить, написав для $u_0$
\begin{displaymath}
\sum_{n=0}^{N}S_nQ(f^{(n)}(u_0) - ku^{(n)}(u_0)) = \sum_{n=0}^N J_n(u_0),
\end{displaymath}
где
\begin{displaymath}
J_n(u_0) = S_nQ(f_n-ku_n)
\end{displaymath}

Корректность и сходимость вышеуказанных допущений и процедур устанавливает теорема.
\newline
\textit{Теорема.} Пусть у диссипативного дифференциального уравнения есть компактный аттрактор $\mathds{A}\subset\mathds{R}^d$. Тогда для любого $\epsilon>0$ 
существуют $k$ и $N$ такие, что
\begin{displaymath}
\left|Qu_* - \sum_{n=0}^N J_n(u_*)\right| \le \epsilon, \:\: \forall u_* \in \mathds{A}
\end{displaymath}

Доказательство теоремы приведено в \cite{1988PhyD...32..163F}.

\section{Пример. Система Лоренца}
В самой статье \cite{1988PhyD...32..163F} приведен пример построения подобных аппроксимирующих множеств для аттратора Лоренца. К сожалению, в нем допущена ошибка при
 вычислении самого первого многочлена. В самом деле, на стр. 176 авторы переписывают уравнение (4.4), которое выглядит как
\begin{displaymath}
\xi = I_0(u), \: \xi=x-y
\end{displaymath}
в виде
\begin{displaymath}
(1+s_0k)\xi+s_0y+s_0(\xi+y)z=0
\end{displaymath}
используя тот факт, что
\begin{displaymath}
I_0(u)=J_0(u)=-s_0(y+k\xi)-s_0(\xi+y)z
\end{displaymath}
и это действительно верно, однако заменяя в переписанном тождестве (4.4) статьи \cite{1988PhyD...32..163F} переменную $\xi$ на $x-y$, авторы получают уравнение
\begin{displaymath}
(1+s_0k)x+(1+s_0(k-1))y+s_0zk=0,
\end{displaymath}
в то время как подстановка на самом деле дает
\begin{displaymath}
(1+s_0k)(x-y)+s_0y+s_0(x-y+y)z=(1+s_0k)x-(1+s_0(k-1))y+s_0xz=0,
\end{displaymath}
а в силу того, что все последующие значения $I_k$ могут использовать полученное выражение, которое, как легко убедиться, 
неверно, нельзя быть полностью уверенным в том, что соответствующие вычисления проведены корректно.
\newline
Поэтому рассмотрим пример с системой Лоренца подробнее.
\newline
Задача состоит в оптимальном выборе $A, k, N, \delta$. Для начала исходная система Лоренца (для классических параметрах $\sigma=10, b=8/3,r=28$), выглядящая как
$$ \left\{
\begin{aligned}
& x' + \sigma x - \sigma y = 0,\\
& y' + y + xz = 0,\\
& z' + bz - xy + br = 0.\\
\end{aligned}
\right. $$
была переписана в виде
$$ \left\{
\begin{aligned}
& \xi' + \sigma \xi - (\xi+y) z = 0,\\
& y' + y + xz = 0,\\
& z' + bz - xy + br = 0.\\
\end{aligned}
\right. $$
с целью сделать матрицу $A$ удобной ($\xi=x-y$).
\begin{displaymath}
A=\begin{pmatrix} \sigma & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & b \end{pmatrix},
\end{displaymath}
в качестве $\Lambda$ взяли $\sigma=10$, соответственно проектор стал проектором на первую координату. Очевидно
\begin{displaymath}
R_0=\begin{pmatrix} 0\\ 0\\ br \end{pmatrix},\:\: R_1(u) = \begin{pmatrix} -y\\ 0\\ 0 \end{pmatrix},\:\: R_2(u) = \begin{pmatrix} -(\xi+y)z\\ (\xi+y)z\\ -(\xi+y)z \end{pmatrix},
\end{displaymath}
Далее авторы строят первые несколько многочленов, аппроксимирующие аттрактор, но сейчас это проще сделать с помощью символьных вычислений, потому что, например, 4-й аппроксимирующий многочлен 
выглядит следующим образом: $2.88016315393237*x - 2.0010612729*y - 2.001715671*z + 0.000021185*z*((61*x)/25 - (21289*y)/1375 - (82167*z)/2750 + (61*x*z)/25 + 61/125) - 0.000014919*(x*z - (57*z)/10 + 1/5)*((61*x)/25 - (21289*y)/1375 - (82167*z)/2750 + (61*x*z)/25 + 61/125) + 0.0000149194427*((61*y)/25 + (61*z)/25)*((57*x*z)/10 - (3249*z)/100 - x*(x*z - (57*z)/10 + 1/5) + (25*z*((61*y)/25 + (61*z)/25))/11 + 57/50) + 0.00010776845*x*z + 0.0000175098*x*(x*z - (57*z)/10 + 1/5) + 0.000004973147591*z*((21289*x)/1375 - (6406586*y)/75625 - (78081769*z)/302500 + (82167*x*z)/2750 - (61*x*(x*z - (57*z)/10 + 1/5))/25 + (61*z*((61*y)/25 + (61*z)/25))/11 + 82167/13750) - 0.0000021881849*x*((185193*z)/1000 + (25*z*((61*x)/25 - (21289*y)/1375 - (82167*z)/2750 + (61*x*z)/25 + 61/125))/11 - (3249*x*z)/100 + (57*x*(x*z - (57*z)/10 + 1/5))/10 + x*((57*x*z)/10 - (3249*z)/100 - x*(x*z - (57*z)/10 + 1/5) + (25*z*((61*y)/25 + (61*z)/25))/11 + 57/50) - (285*z*((61*y)/25 + (61*z)/25))/22 + (50*((61*y)/25 + (61*z)/25)*(x*z - (57*z)/10 + 1/5))/11 - 3249/500) + 0.00000932166784*x*((57*x*z)/10 - (3249*z)/100 - x*(x*z - (57*z)/10 + 1/5) + (25*z*((61*y)/25 + (61*z)/25))/11 + 57/50) - 0.00003979512589*z*((61*y)/25 + (61*z)/25) + 0.00004237121748*((61*y)/25 + (61*z)/25)*(x*z - (57*z)/10 + 1/5) + 0.000021553691049$.
\newline
В рамках данной работы была написана программа на $Matlab$, которая генерировала полиномы для аппроксимации и строила поверхности, на которых принимаются значения близкие к 0.
\newline
Для аттрактора Лоренца первые 7 аппроксимаций выглядят следующим образом:
\newline
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori0.jpg}}
\caption{Аппроксимация 0}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori1.jpg}}
\caption{Аппроксимация 1}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori2.jpg}}
\caption{Аппроксимация 2}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori3.jpg}}
\caption{Аппроксимация 3}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori4.jpg}}
\caption{Аппроксимация 4}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori5.jpg}}
\caption{Аппроксимация #5}
\end{figure}
\begin{figure}[H]
\center{\includegraphics[scale=0.35]{lori6.jpg}}
\caption{Аппроксимация #6}
\end{figure}
Действительно видно, что постепенно аппроксимации становятся сложнее с точки зрения фигур и ближе к аттрактору, образовывая все больше точек пересечения.
\section{Пример. Система Ресслера}
В качестве практического применения данного алгоритма была выбрана так же система Ресслера. Однако 
у этой системы есть только инвариантное множество, которое часто называют аттрактором, но аттрактором как таковым
оно не является в силу того, что в любой окрестности этого множества есть траектории, уходящие на бесконечность. 
Попробуем построить аппроксимацию инвариантного множества системы Ресслера, которая описывается дифференциальным уравнением
$$ \left\{
\begin{aligned}
& x' + y + z = 0,\\
& y' - x - ay = 0,\\
& z' - b - z(x-c) = 0.\\
\end{aligned}
\right. $$
В качестве параметров взяты стандартные $a=b=0.2, c=5.7$. Эту систему пришлось переписывать чуть более сложным способом, получилось более громоздко, нежели в предыдущем примере:
$$ \left\{
\begin{aligned}
& \xi' + 6\xi -7.64\eta-b-z\frac{\xi-\eta-z}{-0.44}-1.74z = 0,\\
& \eta'+\eta-a\frac{\xi-z-1.44\eta}{-0.44}+z=0,\\
& z' - b - z\frac{\xi-\eta-z}{-0.44} + zc=0.\\
\end{aligned}
\right. $$
При этом $\xi=-1.44x+y-z, \eta=y-x$.
Зато 
\begin{displaymath}
A=\begin{pmatrix} 6 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 5.7 \end{pmatrix},
\end{displaymath}
а
\begin{displaymath}
R_0=\begin{pmatrix} -0.2\\ 0\\ -0.2 \end{pmatrix},\:\: R_1(u) = \begin{pmatrix} -7.64\eta-1.74z\\ -ay+z\\ 0 \end{pmatrix},\:\: R_2(u) = \begin{pmatrix} -z\frac{\xi-\eta-z}{-0.44}\\ 0\\ -z\frac{\xi-\eta-z}{-0.44} \end{pmatrix},
\end{displaymath}
Далее с помощью той же программы на $Matlab$, что и из предыдущего пункта были сгенерированы первые 11 полиномов, после чего для каждого посчитали нулевые поверхности.
Начальные данные $k=0, \Lambda=6, \delta_0=10^{-6}$. Получены следующие аппроксимации:

\section{Символьные вычисления}
Как уже было сказано в предыдущих разделах, удобно использовать символьные вычисления для применения данного алгоритма. В самом деле, в силу того, 
что все фигурирующие переменные и множества являются полиномиальными функциями от $u_0$, символьные вычисления позволяют быстро и 
сразу в электронном виде получать искомые полиномы. Рассмотрим программу, считающую первые 3 приближения для системы Лоренца и 
записывающую полученные полиномы в виде функции $Matlab$ для дальнейшего использования.
\lstset{ %
language=Matlab,                 % выбор языка для подсветки (здесь это С)
basicstyle=\small\sffamily, % размер и начертание шрифта для подсветки кода
numbers=left,               % где поставить нумерацию строк (слева\справа)
numberstyle=\tiny,           % размер шрифта для номеров строк
stepnumber=1,                   % размер шага между двумя номерами строк
numbersep=5pt,                % как далеко отстоят номера строк от подсвечиваемого кода
backgroundcolor=\color{white}, % цвет фона подсветки - используем \usepackage{color}
showspaces=false,            % показывать или нет пробелы специальными отступами
showstringspaces=false,      % показывать или нет пробелы в строках
showtabs=false,             % показывать или нет табуляцию в строках
frame=single,              % рисовать рамку вокруг кода
tabsize=2,                 % размер табуляции по умолчанию равен 2 пробелам
captionpos=t,              % позиция заголовка вверху [t] или внизу [b] 
breaklines=true,           % автоматически переносить строки (да\нет)
breakatwhitespace=false, % переносить строки только если есть пробел
escapeinside={\%*}{*)}   % если нужно добавить комментарии в коде
}
%\begin{verbatim}
\begin{lstlisting}[label=some-code2,caption=Generating Polynoms]
% initializing parameters
delta = 0.12;
k = 0;
eigenvalue = 10;

% introducing symbolic variables for u_0
x = sym('x');
y = sym('y');
z = sym('z');

% introducing matrix A and vector R0
% functions R1 and R2 are defined in separated files
% as well as function s_n
A = [10 0 0; 0 1 0; 0 0 8/3];
R0 = [0; 0; 8*28/3];

% introducing symbolic variable for u_i and f_i
u0 = sym('u0');
u1 = sym('u1');
u2 = sym('u2');
f0 = sym('f0');
f1 = sym('f1');
f2 = sym('f2');

% symbolically computing u_i and f_i using x, y and z
u0 = [x - y; y; z];
u1 = -A * u0 - R0 - R1(u0) - R2(u0, u0);
u2 = -A * u1 - R1(u1) - R2(u1, u0) - R2(u0, u1);

f0 = R0 + R1(u0) + R2(u0, u0);
f1 = R1(u1) + R2(u1, u0) + R2(u0, u1);
f2 = R1(u2) + R2(u2, u0) + 2 * R2(u1, u1) + R2(u0, u2);

% calculating s_i
s0 = s_n(delta, 0, k + eigenvalue);
s1 = s_n(delta, 1, k + eigenvalue);
s2 = s_n(delta, 2, k + eigenvalue);

% symbolically computing j_k
j0 = s0 * Q(f0 - k * u0);
j1 = s1 * Q(f1 - k * u1);
j2 = s2 * Q(f2 - k * u2);

% symbolically computing polynoms we are trying to obtain
i0 = j0 - Q(u0);
i1 = j0 + j1 - Q(u0);
i2 = j0 + j1 + j2 - Q(u0);

% writing obtained polynoms into the file created as matlab function
fid = fopen('polynom_i0.m','wt');
i0 = char(i0);
fprintf(fid, 'function answer = polynom_i0(x, y, z)\n answer=');
fprintf(fid,'%s',i0);
fprintf(fid, ';\n end\n');
fclose(fid);

fid = fopen('polynom_i1.m','wt');
i1 = char(i1);
fprintf(fid, 'function answer = polynom_i1(x, y, z)\n answer=');
fprintf(fid,'%s',i1);
fprintf(fid, ';\n end\n');
fclose(fid);

fid = fopen('polynom_i2.m','wt');
i2 = char(i2);
fprintf(fid, 'function answer = polynom_i2(x, y, z)\n answer=');
fprintf(fid,'%s',i2);
fprintf(fid, ';\n end\n');
fclose(fid);

% now we can use this functions for calculating surfaces f(u)=0
\end{lstlisting}
Можно заметить, что в дальнейшем для генерации полиномов более высоких порядков потребуются достаточно большие по порядку 
биномиальные коэффициенты, которые достаточно тяжело выписывать вручную. Для генерации строчек вида 
\begin{displaymath}
u_{k+1}=-Au_k-R_1(u_k)-\sum_{j=0}^k C_k^j R_2(u_{k-j}, u_j)
\end{displaymath}
можно написать отдельный скрипт, автоматизировав таким образом процесс полностью, от исследователя будет требоваться лишь ввод 
начальных параметров системы и констант $k, N ,\delta$.

% У заключения нет номера главы

\section*{Заключение}
Данный метод позволят аппроксимировать аттракторы с хорошей точностью, предоставляя свободу выбра сразу нескольких переменных. Подобные алгебраические множества 
просты для вывода и построения, а потому степень аппроксимации зависит только от доступных вычислительных мощностей.
В то же время столь гибкая структура аппроксимации позволяет вводить дальнейшие структуры для изучения сложных множеств -- например стратификацию.
В данной работе продемонстрирован алгоритм вместе с его работой, которую можно назвать приемлемой. Для более точных расчетов потребовалось бы больше априорных вычислений 
и больше вычислительных мощностей непосредственно при исполнении алгоритма.
\bibliographystyle{ugost2008ls}
\bibliography{AttractorApproximation.texdiploma.bib}
\end{document}